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ABSTRACT 

The SCUBA instrument on the James Clerk Maxwell Telescope has already had an 
impact on cosmology by detecting relatively large numbers of dusty galaxies at high 
redsliift. Apart from identifying well-detected sources, such data can also be mined 
for information about fainter sources and their correlations, as revealed through low 
level fluctuations in SCUBA maps. As a first step in this direction we analyse a small 
SCUBA data-set as if it were obtained from a Cosmic Microwave Background (CMB) 
differencing experiment. This enables us to place limits on CMB anisotropy at 850 /im. 
Expressed as Qa a t, the quadrupole expectation value for a flat power spectrum, the 
limit is 152/iK at 95 per cent confidence, corresponding to C 0 <355 fiK for a Gaus¬ 
sian autocorrelation function, with a coherence angle of about 20-25 arcsec; These 
results could easily be reinterpretted in terms of any other fluctuating sky signal. This 
is currently the best limit for these scales at high frequency, and comparable to limits 
at similar angular scales in the radio. Even with such a modest data-set, it is possible 
to put a constraint on the slope of the SCUBA counts at the faint end, since even 
randomly distributed sources would lead to fluctuations. Future analysis of sky cor¬ 
relations in more extensive data-sets ought to yield detections, and hence additional 
information on source counts and clustering. 

Key words: cosmic microwave background - cosmology: observations - methods: 
data analysis - infrared: galaxies 


1 INTRODUCTION 

Since COBE discovered the existence of anisotropy in 
the Cosmic Microwave Background Radiation (CMB) at 
large angular scales, balloon and ground based experiments 
have detected anisotropy at a range of smaller scales (see 
e.g. Smoot & Scott 1998). All cosmological models predict 
that the very smallest scales should be free of primordial an¬ 
isotropy, because photons in small-scale overdensities which 
entered the horizon early have been able to random walk 
out of the potential wells, and also because fluctuations on 
scales below the thickness of the last scattering surface are 
suppressed. 

Nevertheless, secondary anisotropies can be generated 
through a wide variety of physical processes occurring at 
redshifts < 1000. For an overview of some of these effects 
see Bond (1996). There is new motivation from the detec¬ 
tion of the Far Infrared Background (FIB, Puget et al. 1996, 
Schlegel, Finkbeiner & Davis 1998, Fixsen et al. 1998, 
Hauser et al. 1998), as well as recent SCUBA results, that 
dusty galaxies at high redshifts could be of greater signif¬ 
icance than previously assumed - possibly the dominant 


source of CMB anisotropy at arcsecond scales, and certainly 
the dominant source of sky fluctuations at ~ 1 mm, at least 
out of the galactic plane. The smooth FIB is expected to 
break-up into sources, plus fluctuations due to unresolved 
sources, and correlations between sources on the sky. Weak 
Sunyaev-Zel’dovich increments, caused by hot gas in viri- 
alised structures, may also contribute to anisotropies in the 
sub-millimetre, and of course other physical effects could 
also be important. We attempt here to use SCUBA data to 
place limits on general fluctuations within the framework of 
CMB anisotropies, i.e. as intensity fluctuations on a 2.728 K 
blackbody spectrum; but everything can be reinterpretted 
in terms of fluctuations of some other spectral component. 

Experiments probing similar angular scales have al¬ 
ready placed upper limits of anisotropy at a variety of wave¬ 
lengths, and a summary can be found in Table 1. The results 
are generally quoted in terms of the most sensitive GACF 
(which will be described in Section 4), although a few ex¬ 
periments simply give a limit to the variance of the tem¬ 
perature measurements. Limits set at radio wavelengths are 
markedly more sensitive than those in the millimetre region. 
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Figure 1 . The arrangement of the 37 bolometers in the long 
wavelength array. The squares of 9 points are the jiggle positions 
that create an effective beam of 14.7 arcsec, which are shown as 
the larger circles on the main plot. The solid jiggle pattern shows 
the bolometer positions for the first integration. Due to earth ro¬ 
tation, the entire array rotates throughout the course of the run, 
and we have plotted the position of the bolometers for the last 
integration using open circles. The inset illustrates the three posi¬ 
tions for the three beam chop that constitutes a single integration. 
The entire array is chopped in azimuth with this double-difference 
pattern. 

However, the expected foreground signals at low frequency 
are expected to be entirely different to those which are likely 
to dominate fluctuations in the sub-millimetre sky (namely 
dusty galaxies and hot cluster gas). The most relevant pre¬ 
vious measurement is the earlier .JCMT limit of Church, 
Lasenby & Hills (1993), using the single bolometer which 
was the predecessor to SCUBA. 


2 SCUBA OBSERVATIONS 

The observations were conducted with the Submillimeter 
Common-User Bolometer Array (Cunningham et al. 1994; 
Gear & Cunningham 1995; Lightfoot et al. 1995) on the 
James Clerk Maxwell Telescope. SCUBA contains a number 
of detectors and detector arrays cooled to 0.1 K that cover 
the atmospheric windows from 350 pm to 2000 pm. 

The arrays have a hexagonal arrangement of pixels, with 
feeds about two beamwidths apart in the focal plane. The 
two array detectors provide an instantaneous field-of-view 
of 2.3 aremin and can be used simultaneously. They consist 
of the 91 element Short-wave array, which we used with the 
450 pm filter, and the 37 element Long-wave array, which 
we used at 850 pm. An illustration of the observing strat¬ 
egy is given in Fig. 1, and described in more detail below. 
Essentially it is a double difference experiment, with some 
additional complications. 


During an integration, the telescope is ‘jiggled’ in a 3 
x 3 square pattern (with 2 arcsec offsets, shown in Fig. 1). 
This improves the photometric accuracy, reducing the im¬ 
pact of pointing errors by averaging the source signal over 
an area slightly larger than the beam. The effective half¬ 
power beamwidths for the 450 and 850 pm pixels are then 
7.5 arcsec and 14.7 arcsec respectively. 

Whilst jiggling, the secondary was chopped at 7.8125 Hz 
between two positions, A and B, separated here by 90 arcsec 
in azimuth, yielding the single difference measurement Ib — 

1 a- This is followed by a second single difference, Ic~Ib , af¬ 
ter the nod to the other side of B. By differencing these two 
differences, we obtain a series of double differences, with the 
signal defined by S = —Ia + “21b — Ic- For some of the sys¬ 
tematic checks we will consider the data on a jiggle by jiggle 
basis. However, the data are generally treated by lumping 
the 9 jiggles together, yielding an 18 second double differ¬ 
ence for each bolometer, which we refer to as a single ‘inte¬ 
gration’. Considered in this way the set-up is the traditional 
triple beam arrangement common in CMB experiments. 

SCUBA can also use a 64 point jiggle in order to fill in 
the blank space between the pixels (‘mapping mode’), but 
for this initial study, we use ‘photometry’ observations at 450 
and 850 pm. Although this provides an undersampled map 
at both wavelengths, it is more straightforward to analyze. 

On 1997 December 3, five ‘scans’ of 900 seconds, con¬ 
sisting of 50 integrations each were taken. The central pixel 
of SCUBA was fixed on a point source in an otherwise blank 
field, for which we were interested in obtaining sub-mm pho¬ 
tometry (the lensed AGN B1933+503, Chapman, Fahlman 
& Scott, in preparation). The other pixels rotated relative to 
the sky through the period of the integration; we show the 
positions of the first and last integrations in Fig. 1. Pointing 
was checked hourly on the blazar 2036-419 and a sky-dip 
was performed between each 15 minute scan to measure the 
atmospheric opacity. The rms pointing errors were below 

2 arcsec, while the average atmospheric zenith opacities at 
450 pm and 850 pm were fairly stable with r being 0.51 and 
0.12 respectively. However, there were some short time-scale 
variations, presumably due to water vapour pockets blowing 
over at high altitude, which caused some parts of the data¬ 
set to be noisier (see Fig. 2. The observations were largely 
reduced using the Starlink package SURF (Scuba User Re¬ 
duction Facility, Jenness & Lightfoot 1998). Spikes were first 
carefully rejected from the double difference data. The data 
were then corrected for atmospheric opacity and calibrated 
against Saturn and the compact HII region K3-50, which 
were also observed during the same observing shift. The 
850 pm calibrations agreed with each other and also the 
standard gains to within 10 per cent. However, at 450 pm, 
K3-50 is extended and variable, and is not a good calibra¬ 
tion source, while the Saturn 450 pm calibration agreed with 
standard gains to within 25 per cent. 


3 DATA REDUCTION 

It is necessary to identify and remove any non-astronomical 
signals from the data. One of SCUBA’s many strengths is its 
ability to redundantly measure the strongest of these con¬ 
taminants: atmospheric emission. To a lesser extent, cosmic 
ray hits influence the data, but these are readily identified 
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Table 1. Summary of previous small angular scale CMB limits. 


A 

Angular Scale 

A T/T 

Instrument 

Reference 

800 pm 

17 arcsec 

< 146 x 

10~ 5 

JCMT 

Church et al. (1993) 

1250 pm 

30 arcsec 

< 14 x 

10 -5 

SEST & IRAM 

Andreani (1994) 


70 arcsec 

< 24 x 

10~ 5 




140 arcsec 

< 19 x 

10' 5 



1300 pm 

12 arcsec 

< 18 x 

10~ 5 

IRAM 

Kreysa & Chini (1989) 

2110 pm 

66 arcsec 

< 2.1 x 

10" 5 

SuZIE 

Church et al. (1997) 

3400 pm 

10 arcsec 

< 9 x 

10~ 5 

IRAM 

Radford (1993) 

15,000 pm 

156 arcsec 

< 1.7 x 

10~ 5 

OVRO 

Readhead et al. (1989) 

19,700 pm 

20 arcsec 

< 2.3 x 

10~ 5 

Ryle 

Jones (1997) 

35,000 pm 

30 arcsec 

< 2.3 x 

10~ 5 

ATCA 

Subrahmanyan et al. (1998) 


60 arcsec 

< 1.6 x 

10" 5 




120 arcsec 

< 2.5 x 

10” 5 



35,000 pm 

6 arcsec 

12.8 x 

10~ 5 

VLA 

Partridge et al. (1997) 


18 arcsec 

< 4.8 x 

10" 5 




60 arcsec 

< 2.0 x 

10” 5 




by the spikes they leave in the timestream, and are there¬ 
fore easily removed from the data (more details in the next 
section). Finally, we test for other correlations in the data 
that indicate a common signal, such as cross-talk between 
bolometers. 

In order to discuss our analysis procedure, we adopt the 
following notation. The indices b, s, i, and j denote bolome¬ 
ter, scan, integration, and jiggle number respectively, each 
starting at 1. The variable N subscripted with one of these 
indices represents the total number of the quantity. For 
this particular data-set they are: Nb = 37 or 91 (for the 
450 pm and 850 pm channels respectively); N s = 5; Ni = 50; 
and Nj = 9. It will also be useful to introduce the variable 
k = (s— 1 )Ni + i, which simply indexes a particular integra¬ 
tion. 


3.1 Removing the atmospheric signal 

We plot the signal timestream for two representative 
bolometers from each of the short and long-wavelength ar¬ 
ray in Fig. 2. It is clear that the output is highly correlated 
in time between bolometers, even at different wavelengths. 
Furthermore, a detailed inspection of all the timestreams 
shows that the correlation is strong even for bolometers at 
opposite ends of the array, indicating that a common at¬ 
mospheric signal subtends an angle greater than 2 arcmin. 
Order of magnitude considerations suggest that the size of a 
relevant patch is perhaps 1000 arcsec (e.g. Jenness, Lightfoot 
& Holland 1998). Much of this signal is removed in the pro¬ 
cess of chopping and nodding, but some atmospheric noise 
inevitably remains (e.g. Duncan et al. 1995). It is reasonable 
to use the correlation across the array to calculate a com¬ 
mon atmospheric signal that can be removed from the data. 
We will now consider two separate methods of removing the 
atmospheric signal from the long wavelength data, namely 
use of the average across the long wavelength array, or of 
the independent average from the short wavelength array. 



Figure 2. The raw timestreams from 2 bolometers on each of 
the 450 //in and 850 //rn arrays. The central pixels at 450 //rn and 
850 pm are shown in panels b) and d). The pixels plotted in panels 
a) and c) were chosen to be far from the central pixel to help 
illustrate the correlated signal caused by the atmosphere. Vertical 
dashed lines separate the 5 ‘scans’ (sets of 50 integrations). 


3.1.1 Using the 850pm data to remove atmospheric 
signals 


The raw data contains the 2 second double difference signal, 
Vbkj which we use to compute a mean sky signal on a jiggle 
by jiggle basis: 


Mkj = 


1 

N b - N b E 


N b 

^ ] v bkj■ 
b=l,b<£E 


(1) 


Here N b is the number of bolometers in the set E, which 
is a list of all bolometers that are excluded from this sky 
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calculation. These include the central pixel (which is ‘con¬ 
taminated’ by B1933+503) and any bolometers that ex¬ 
hibit excessive noise. These were chosen by looking at the 
sky-corrected, integrated signal where the mean was com¬ 
puted leaving out only the central pixel. Those excluded 
show an integrated variance about twice as high as the 
other bolometers in their channel. For the 850 pm data, 
7 bolometers were removed (E= {12,19, 22, 23, 24, 32, 37}) 
and the integrated signal changed by only 1 pV, with 
a comparable reduction in the error bar. At 450 /im, 
E= {46, 87,8, 59, 62, 64, 65, 66, 67, 69, 70, 71, 78, 79} and the 
integrated signal changed by no more than 4 pV. The sky- 
corrected signal is then simply 

Vbfcj — Vbkj Mfcj . (2) 


The 9 jiggles at each integration are now binned to¬ 
gether by taking a direct mean to form the signal Vbk ■ Be¬ 
cause the readout noise is very stable for a given scan, we 
can assume that the noise in each integration is drawn from 
the same random distribution. If we further assume that this 
noise is uncorrelated from integration to integration, we can 
compute the mean and variance of the sky-corrected signal 
for each scan and each bolometer using the data, via the 
scatter of the points within a scan: 


V 6s 



N i 


^ ^ VbsiQbsi 

i= 1 


( 3 ) 


and (aY s f = 


1 


Ni-Nf - 


— (Vf , s i — Vbs) 


Qbs 


( 4 ) 


Integrations k = 1-40 and k = 101 140 (the first parts of 
scans 1 and 3) were taken while clouds were obscuring the 
target region (see Fig. 2), rendering them unusable because 
the rapidly changing opacity cannot be easily characterised. 
To account for these data, we introduce the quantity N£, 
which is the total number of bad integrations in a particular 
scan for a given bolometer, and the quantity qw, which 
takes the values 1 and 0 for good and bad data respec¬ 
tively. For our particular data set, there are then 80 x Nb = 
2960(7280) integrations at 850(450) /im that are not used in 
the subsequent analysis. 

Anomalous signals (e.g. cosmic ray hits) in the cor¬ 
rected data were removed on a bolometer by bolometer, scan 
by scan basis. Any integration that deviated by more than 
3 <j^ was removed and the variance recalculated. This pro¬ 
cedure was repeated once more to ensure that any statis¬ 
tically significant anomalies shadowed by even larger ones 
were removed during the second pass. The total number of 
integrations removed was 31/2 (44/4) for the two passes at 
850(450) fjm. A third pass failed to remove any additional 
points. The removal of noisy sky sections and anomalous 
signals left 68 per cent of the data to be used in subsequent 
analysis. Finally, Fig. 3 plots the integrated signal for each 
bolometer, which is calculated using the weighted mean 


h 


N * Er=i «r 2 


( 5 ) 


Note that the likelihood analysis in Section 4 does not simply 
use this binned version of the data, but allows for full spatial 
correlations on the sky. 



Bolometer 


Figure 3. The integrated signal measured by each bolometer on 
the 450 pm array is shown in the top panel and the 850 pm array 
in the centre panel. The detection in the central pixels, bolometer 
19 at 850 pm and 46 at 450 pm, is readily apparent and due to the 
known source at that position. The variance of each bolometer is 
also plotted separately in the bottom panel, to clearly identify 
noisy bolometers. 


The method of removing atmospheric signal described 
above is very successful, and some variant of it is the usual 
method adopted for reducing SCUBA data. One disadvan¬ 
tage however, is that its use requires knowing in advance 
which pixels contain signal, or using some iterative process 
to remove the sky when low-level signals are present. An¬ 
other disadvantage, peculiar for our purposes here, is that 
this method correlates all of the pixels, since the mean has 
been subtracted from each one. It is possible to take this into 
account in a full likelihood analysis of the fluctuations by es¬ 
sentially also removing a mean from the theory. However, we 
simply ignored these correlations, realizing that they will 
be negligible for the bin sizes used in our analysis. Each 
binned data point uses ~ NbNi integrations in determining 
the mean signal, introducing a tiny correlation that is the 
inverse of this quantity. Also, any correlations in the data 
will show up as a signal in the likelihood plots of Section 4. 
Since we ultimately find that the signal is consistent with 
zero, we can safely assume that our mean subtraction ap¬ 
proach is sufficient. Were this not the case, we would have 
to make the analysis insensitive to the mean using a ma¬ 
trix rotation method (Bunn et al. 1994), or marginalization 
(Bond, Jaffe & Knox 1998). 


3.1.2 Using the 450 pm array as an atmospheric monitor 

We can avoid the correlation problem altogether by using 
an independent estimate for removing the baseline from the 
data. The mean sky signals calculated from equation (1) are 
extremely correlated for the 450 pm and 850 pm data, and 
therefore it is feasible to attempt to subtract the sky us- 
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Figure 4. Using the 450 pm array to remove the atmosphere from 
the 850 pni array. The top panel shows the mean sky signals for 
both arrays over a subset of the timestream, with the 850 firn data 
solid and the 450 fim data dashed. The signals are clearly corre¬ 
lated. The bottom panel demonstrates the small residual between 
the 850 fin r mean and scaled 450 fin l mean. 

ing the independent data from the other channel. This may 
be particularly useful for future SCUBA cosmology studies, 
since for ‘blank’ fields the 450 pm data will generally con¬ 
tain no signal (while the 850 fi m data may contain a contri¬ 
bution from extragalactic sources). Hence this method may 
have general utility for helping to look for low levels of fluc¬ 
tuations in long integrations at 850 fin i, where it becomes 
important not to remove any of the signal along with the 
sky. 

The 450 pm channel is more susceptible to changes in 
opacity, so we divide each scan into 10 sections with 5 inte¬ 
grations each, and perform a least-squares fit of the form 

x 2 = EEK-«+»)] 2 ^ (6) 

k j 

where k runs over the 10 integrations to be fit (excluding, of 
course those integrations dominated by noisy atmospheric 
signal) and j over each jiggle point in the integration. The 
value of c is typically near 0.2 (see also Jenness et al. 1998), 
and does not vary by more than 10 per cent within a scan. 
The offset o is on the order of a few pV, and we found it 
necessary to include it; If we redo the fit and fixo = 0.0, the 
integrated signal tends to be biased lower by approximately 
4 pV. In Fig. 4 we plot the 450 pm and 850 pm mean sky for 
a section of scan 4. Also in the figure we plot the residual be¬ 
tween the 850 pm mean sky and the scaled 450 pm data. To 
appreciate how small this residual is, we take as an example 
the rms of the residual for the 450 jiggles in scan 4, which 
turns out to be 100 pV. The rms of the sky-corrected signal 
for a typical bolometer in scan 4 using either the 850 pm or 
450 pm mean is about 380 pV. 

After forming the new mean signal, we subtract it from 


Figure 5. The difference in the integrated signals between 850 pm 
data sky-subtracted using the 850 and 450 pm array averages. 

the 850 pm data and calculate the integrated signal. As ev¬ 
ident in Fig. 5, the signal level changes by at most 4 pV 
compared to the results obtained by using the 850 pm data 
to subtract a sky signal (top panel). There does not appear 
to be any bias introduced due to this method, although 
the variance is systematically higher by about 0.5 pV or 
0.12mjy (bottom panel). The increased variance is consis¬ 
tent with a signal with an additional independent noise term 
about one-third the size of the main noise term, as we can 
predict from the residual of the 850 and 450 pm mean signal 
level. 

Therefore this method invariably introduces more un¬ 
certainty in the integrated signal, but may be the best op¬ 
tion when there is, for example, extended structure in the 
850 pm map, and nothing but noise in the 450 pm. In this 
case, however, care must be taken to avoid removing a DC 
level in the 850 pm data via the parameter o, which is a sys¬ 
tematic effect that could be calibrated by analyzing other 
data sets that are largely free of sources in both channels. 
For our data set, this additional variance is only 5 per cent, 
and influences our likelihood fits in Section 4 by roughly this 
amount. 

3.1.3 Testing for a plane through the bolometer array 

It is plausible that the two dimensional bolometer array ex¬ 
hibits evidence of a systematic temperature gradient, caused 
either by atmospheric signal or by sources even more local 
to the telescope. This might introduce extra variance into 
the data, which could erroneously be identified with larger 
scale signal. To clarify this, we attempt to fit the data for 
each array using 

P(x, y) = A + Bx + Cy (7) 

where x and y are the differences in altitude and azimuth 
respectively from the central beam. 
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Having found the best fitting plane at each time in¬ 
terval, we want to assess whether it is real (i.e. due to the 
atmosphere or some other systematic effect) or just statis¬ 
tical (i.e. due to stochastic noise). The simplest approach 
is to test the hypothesis that a for a given integration there 
exists a correlation between the best fit plane in the 450 and 
850 /tm data, as would be expected if the plane is a genuine 
atmospheric effect. To implement this, we apply a Spearman 
rank correlation test to the following two statistics: 
the normalised height of the edge of array with y=0, 


A'(A) 


B x 70" 
A 


and at x=0, 


y(A) 


C x 70" 
A 


( 8 ) 

(9) 


Here 70 arcsec is the approximate radius of the SCUBA 
array. 

The results are somewhat more complicated than we ex¬ 
pected. Performing the analysis for each integration, we find 
a highly statistically significant correlation of r = 0.44 be¬ 
tween planes in the azimuthal direction, but with a relatively 
small slope: (F(850/um) = 0.01; Y(450/um) = 0.04). This 
small residual slope is not surprising, considering that a dou¬ 
ble difference measurement inherently removes gradients, 
but the significance of the correlation (formally P — 2 x 10~ 9 
was unexpected). 

The correlation and significance parameters of the test 
performed on the altitude direction are comparable (r = 
0.33, with a probability of P = 1 x 10 —5 ), however, the 
slopes of the planes are, on average, more pronounced: 
X(850/iin) = -0.11; X(450/xm) = 0.01. 

These results certainly indicate that there is a residual 
slope across the data, most likely as a result of atmospheric 
variation. The fit for the 450 /im plane is generally better (as 
verified by inspecting of the reduced chi-square statistic), 
because there are more bolometers in the short-wavelengtli 
array. Particularly at 850 /tm the best-fitting plane for a par¬ 
ticular integration will be largely a stochastic variation, al¬ 
though taking off these planes will statistically remove some 
atmospheric signal. 

We also looked at whether there were correlations in 
the planes when we binned integrations together, in order to 
bring down the stochastic variation in the planes. We found 
that the highest correlation we could obtain was r = 0.60 
when the bin size was 5 integrations. However, there are then 
only 170/5 = 34 planes to use in the rank correlation test, 
making the probability derived from the rank correlation 
questionable. Indeed, this fit did a poor job of cleaning up 
the signal, because the rms after sky subtraction using this 
plane was twice as large as the single integration plane. With 
the data binned into 5 integrations we found the azimuthal 
correlation between the 450 and 850 /mi planes to be not 
very significant, while the altitude slope correlation was still 
very highly significant. 

The conclusion is that some evidence of residual sky- 
planes exists in our data set, perhaps stronger in the altitude 
direction than the azimuth direction. These effects, however, 
are quite small, and are almost negligible for our analysis of 
fluctuations. Nevertheless we will regard the data-set with 
sky-planes removed as providing our best estimate of any 


residual fluctuations. Clearly this issue merits further inves¬ 
tigation with more comprehensive SCUBA data-sets. 


3.2 Testing correlations between the bolometers 

Systematic effects might introduce correlations between 
bolometers in the time domain, while genuine astronomical 
signal shows up as correlations in the space domain. There 
are three types of potential correlations introduced by the in¬ 
strument for which we would like to test: cross-talk between 
the bolometers grouped 16 per A/D card, crosstalk between 
exceptionally noisy bolometers and their neighbours (listed 
in Section 3.1.1), and possibly spill-over from the central 
pixel which contains the signal from B1933+503. 

The cross-correlation coefficient, rxY of the sky- 
corrected signal was taken between each pair of data points. 
Here, we are using the 170 18-second integrations calculated 
previously. We find that |rxy | < 0.3 except for bolometers 
32 and 36, which have a fairly strong positive correlation of 
0.5. In fact these two bolometers are adjacent to each other 
in the array and on the same A/D card. In any case, we 
find that removing data from this pair (or any other pair) of 
bolometers affects the final results by no more than 20 per 
cent. However, we argue that the correlation is not strong 
enough to say with certainty that cross-talk is the cause. It 
is of course possible that both bolometers are observing an 
unresolved source with a flux slightly below the noise level 
(we will return to this issue in Section 5, when we discuss 
the impact of unresolved point sources on the CMB fluctu¬ 
ations). For that reason we do not exclude these bolometers 
from our analysis. 


3.3 Calibration and conversion of data into 
Thermodynamic temperature 


Calibration was performed using sources of known flux dur¬ 
ing the same observing period. Specifically we used Saturn 
and the HII region K3-50. Corrections were also made for 
the zenith opacity and the elevation angle. As discussed in 
Section 2, these calibrations agreed closely with the ‘stan¬ 
dard gain’ for the system at 850/tm. At 450/tm the agree¬ 
ment was not so good, but the shorter wavelength data were 
less important for our analysis in any case. We list in Table 2 
the conversions between voltage and flux density which we 
used. 

The measured intensity of a blackbody at a thermody¬ 
namic temperature, Tcmb, is given by 


h = 


2 hu 3 1 

q 2 gx — ^ 


( 10 ) 


where x = hv/kTcMB■ To get the flux density in the beam, 
we multiply this by the solid angle of the beam: S v = 
Assuming that the beam is Gaussian, the solid angle is 
= (2n/8 In 2) FWHM 2 ~ 1.133 FWHM 2 . Knowing the in¬ 
tensity at a given frequency v, we can calculate the tem¬ 
perature of the source relative to the CMB by taking the 
derivative of equation (10): 

ATcmb - A ^ CMB ^^-Tcmb, (11) 

which can also be written as 
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ATqmb = 


2 c 2 h 2 sinh 2 (x/ 2 ) 

7„3 t 2 7* 

^ ± CMB x 


AS'cmb- 


( 12 ) 


For full accuracy equation (11) should be averaged over the 
bandpass of the filter. Carrying this out for the 850 pm filter 
(see Holland et al. 1998b) we find the effective frequency 
to be 348.4GHz (or 860.5pm). This corresponds to x = 
6.13 for the 850 pm filter, while the value is 11.9 at 450 pm, 
which is therefore too far out in the Wien tail to have any 
sensitivity to CMB fluctuations. The resulting conversion 
factors between Jansky and Kelvin, for the 450 and 850 pm 
data are given in Table 2. 


4 PLACING CONSTRAINTS ON AN 

UNDERLYING ASTRONOMICAL SIGNAL 


4-2.1 The flat power spectrum (Qnut) 

CMB anisotropies can be described by an expansion of 
spherical harmonics on the sky. The angular power spec¬ 
trum of the amplitude of the spherical modes is given by 
Ci, and completely describes the temperature correlations 
between two spots on the sky: 

C{9ij) = (AT(m)AT(nj)) 

= h £,= a (2* + 1 )C e Pi(hi ■ rij) e-^ +1) " 2 , (13) 

where Pi are the Legendre polynomials and cos dij = iii ■ rij. 
The exponential term at the end accounts for the SCUBA 
Gaussian beam shape, with Gaussian width given by a. 

The correlation model is encoded in the Ct, and is the 
quantity we are trying to fit. For the specific case of a scale- 
invariant Sachs-Wolfe, or ‘flat’ power spectrum, this takes 
the form of 


4.1 Binning the data 

Because of sky rotation, the data for each bolometer can¬ 
not be simply averaged across the entire observing run. In¬ 
stead, the timestream is divided into bins that are some 
reasonable fraction of the beamsize. To choose the bin-size, 
let us restrict our attention to the bolometers on the outer 
ring of the array, which will rotate the most throughout 
the course of the observations. Using the positional infor¬ 
mation in the data, we calculate that the central beam in 
the three-beam measurement moved by roughly 6 arcsec (a 
third of a beamwidth) over the entire run. However, the two 
‘off’ beams moved this same distance in just one scan. In 
addition there is a delay between successive scans, so bin¬ 
ning more than one scans’ worth of data would result in a 
smeared beam. Binning the data any finer than one scan 
would simply result in a more unwieldy and noisier data¬ 
set, with no additional information. This would not be the 
case if the bolometer noise varied throughout a scan, but 
for our data-set constant noise in each bolometer is a good 
approximation. 

With these considerations, plus removing those sections 
of data contaminated by atmospheric emission, we obtain 
175 binned spatial points to use in the likelihood fits (5 
spatial points for each of the 35 usable bolometers). Since 
there is potentially significant information contained in the 
off-diagonal correlations between sky positions, we wanted 
to avoid binning any more coarsely than this. In addition we 
also want to retain the possibility of negative correlations 
introduced between ‘on’ beams and ‘off’ beams when we 
perform our likelihood analysis, as described in the next 
section. 


4.2 Likelihood analysis 

Assuming that the noise and underlying astronomical signal 
are Gaussian distributed, all the information is contained in 
the two-point correlation function, (AT(hi)AT(hj)), where 
h is a unit vector on the sky, and i, j correspond to two data 
points. The traditional approach is to take a parameterized 
theory that describes the signal and perform a Bayesian like¬ 
lihood analysis on the data in order to determine the value 
of the parameters. We consider here the two simplest mod¬ 
els for underlying sky fluctuations, the flat power spectrum, 
and the Gaussian Autocorrelation Function (GACF). 


= 24-7T Qj at 

5 e(i + i)’ 


(14) 


where Qfl at is the amplitude parameter to fit using the data 
(and explicitly is the extrapolation to the quadrupole for a 
flat spectrum). 

For an experiment with N data points, we compute 
equation (13) for each pair of points and construct Cfj, the 
theoretical correlation matrix. It completely describes the re¬ 
lationships between the data, and depends on the model of 
the anisotropy (the Cp), the positions on the sky that the ex¬ 
periment studies (the Legendre polynomials), and the beam 
response function (the exponential). 

These last two terms can be considered as a weighting 
function of the power spectrum, which is called the ‘window 
function’ (White and Srednicki 1994) of the experiment, 

WiiOij) = Pe(rii ■ rij)e- l{l+1)a2 , (15) 


and can be used to judge the range of scales to which an 
experiment is sensitive. 

This is the simplified case for an experiment that mea¬ 
sures the temperature of single spots on the sky. A SCUBA 
data point actually consists of three beams with weights 
tu = 1-1, 2,-1]: 

3 

AT/ = ^w a AT(n“), (16) 

oc=1 

which changes the window function and correlation matrix 
to 

3 3 

WiiOij) = EE- w 0 P e (hC^)e-^ +1 ^ 2 , (17) 

a /3 


3 3 

CS‘ = EE^( C “- 1 W-^))- ( 18 ) 

a /3 


For estimating the scale to which the experiment is sen¬ 
sitive, the off-axis contributions (i j) are ignored. Then 
equation (17) reduces to the window function at zero lag: 

W e = e -<v+1)CT2 [3 - 4 P;(cos 7 ) + Pfl cos2 7 )], (19) 

where 7 is the chop amplitude. Fig. 6 is a plot of the zero lag 
window functions relevant for a SCUBA chop of 90 arcsec. 
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Table 2. Conversion factors between thermodynamic temperature, flux density and voltage. 


Wavelength ( fim) 

Convert from 

to 

Multiply by 

Comment 

450 

Voltage 

Flux density 

800 Jy/V 


450 

Flux density 

Temperature 

49,800 AtK/mJy 

FWHM=7.5 arcsec 

850 

Voltage 

Flux density 

240 Jy/V 


850 

Flux density 

Temperature 

568 /rK/mJy 

FWHM=14.7 arcsec 



Figure 6. The zero lag window functions for the double-difference 
SCUBA set-up with a beam throw of 90 arcsec. The high l be¬ 
haviour is set by the beamsize, while the low t shape would be 
different for other beam throws. 


The effective center, (£) is derived by taking the window 
function weighted with a flat power spectrum. 

Once we have the correlation matrix, we can fit for the 
the parameter Qfl at using 

L{Qa at) oc A- exp (d • K _1 (Q fla t) • d T ) , (20) 

With /STy(<2flat) = Cy(Qflat) + Mij (21) 

(Bunn et al. 1994; Srednicki et al. 1994). Here d is the 1 x 
N vector of measured data, and A f is the noise correlation 
matrix. We assume that the data have uncorrelated noise 
signals after sky subtraction, and thus Afij = cno'jSij. 

The matrix inversion can be solved using standard tech¬ 
niques such as singular value decomposition, but it is more 
computationally efficient to use advanced methods, of which 
signal to noise eigenmode decomposition is the most com¬ 
mon (see e.g. Bond 1995; Bunn 1995; Knox 1997), when the 
number of data points is on the order of 200 or so. 

Note that in previous CMB experiments at these scales, 
it was possible to neglect off-axis contributions because of 
the iarge anguiar separation between data points. This is not 
the case for SCUBA, as the pixei separation on the array is 


comparable to the beamsize and is only a factor of 4 smaller 
than the chop amplitude for this data-set. 

A plot of the likelihood function over a range of Qfl at is 
given in Fig. 7. The results indicate that the level of aniso¬ 
tropy is consistent with zero, with an upper limit of 143 /rK 
(95 per cent confidence). This result was obtained using the 
850 fi m mean for sky subtraction, and is shown by the dot¬ 
ted line in the figure. The upper limit when the 450 fim 
mean is used is 164 fiK, approximately 7 per cent higher, 
as expected from the previous discussion, and is shown by 
the short-dashed line. The long-dashed curve in Fig. 7 was 
generated by removing the off-axis contributions to the cor¬ 
relation matrix. This would have given a substantially dif¬ 
ferent (and incorrect) answer, verifying the importance of 
retaining these terms in the analysis. Note that we have in¬ 
tegrated the likelihood using a uniform prior distribution in 
Qflat- Different choices of prior would affect the results some¬ 
what, although the likelihood falls off so rapidly at high Qa a t 
that it is hard to increase the upper limit dramatically by 
any reasonable choice of prior. 

As a test we removed bolometers 32 and 36 from the 
data and redid the analysis. As mentioned earlier, these 
bolometers show stronger evidence of correlation than other 
pairs. The upper limit falls, not surprisingly, to 136 /jK when 
these signals are removed. We also tested the effect of remov¬ 
ing the noisy bolometers, and found that they do not affect 
the Qfl a t limit appreciably. Furthermore, the effect of remov¬ 
ing a plane from the data, as described in Section 3 has only 
a small effect, as we expected from the small gradients that 
we found. Nevertheless, we consider the data-set with the 
sky-planes subtracted from the timestream to yield our best 
estimate of the likehood function, and so we show this curve 
with a solid line in Fig. 7. The upper limit is 152 fiK (95 
per cent confidence). 


4.2.2 Gaussian Auto-Correlation Function (GACF) 

Although largely replaced by the flat power spectrum in 
most recent analyses of CMB anisotropy, it is still useful to 
calculate the result obtained assuming a GACF as a model 
for the correlations. This is particularly true since at these 
angular scales, we do not expect any primary CMB fluctua¬ 
tions, and so any actual signals present could be of any form, 
even one with a preferred correlation scale. For a GACF as¬ 
sumption, the theoretical correlation between two beams of 
Gaussian width a separated on the sky by an angle dij is 
given by 


c(dij) 


c 0 e 2 c [ -e% 

2<J 2 + da 6XP 2(2a 2 + 0 2 ) 


( 22 ) 
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Figure 7. The likelihood function for the flat power spectrum 
model. The solid line shows our best estimate of the likelihood 
function, with subtraction of sky-planes using the 850 fun data, 
and including the full correlation matrix. The dotted line is the 
result obtained using the 850 ftm mean for sky subtraction, while 
the short-dashed line shows the limit using the 450 (im subtrac¬ 
tion. The long dashed line is the 850 fun corrected data had we 
ignored off-axis contributions to the correlation matrix. 


Here Co is the amplitude of the fluctuations, and 9 C is the 
sky coherence angle. 

Again, we construct the correlation matrix for a three 
beam experiment using equation (18), and redo the likeli¬ 
hood analysis. In this case, we need to evaluate the likeli¬ 
hood function over a two dimensional parameter space in 
order to fit for Co and 9 C . 

The result, shown in Fig. 8 for the data sky-subtracted 
using the 850 /./m mean, is C^ 2 < 355 fj,K (95 per cent confi¬ 
dence), at a most sensitive coherence angle of 23 arcsec. To 
be explicit, we integrated using a uniform prior probability 
distribution in the quantity Cq . It is customary to use Cq' 72 
as a measure of the rms temperature sensitivity to compare 
with other experiments; We obtain A T/T < 1.1 x 10~ 4 , 
which is a full order of magnitude lower than the Church 
et al. (1993) result using the precursor to SCUBA, UKT14, 
where the wavelength and coherence angle are are essentially 
identical. Also note that Qflat/Cg' 72 ~ 0.4, and 1 /6 C ~ {£), 
which are the expected comparisons between these analysis 
methods (White and Scott 1994). Fig. 8 makes it clear that 
the data are not strongly correlated, as the confidence limit 
in the GACF is fairly broad in coherence angle. Of course, 
we expect this result based on our Qflat analysis, where the 
most likely correlation was zero. 


Figure 8. The likelihood contours (90, 50, 20 and 10 per cent of 
the peak) for the two parameter GACF model are shown as solid 
lines. The dotted line shows the 95 per cent two dimensional con¬ 
fidence region. The most sensitive coherence length is 23 arcsec, 
where C 0 ' < 355 //K (marginalised 95 per cent confidence limit). 


5 CONTRIBUTION OF UNRESOLVED POINT 
SOURCES TO THE CMB POWER 
SPECTRUM 

In the usual notation, Ci denotes the power spectrum of 
CMB anisotropies on the sky, which is the spherical har¬ 
monic analogue of the Fourier power spectrum on a flat 
plane. Just as in the Fourier case, a population of Poisson- 
distributed point sources will contribute equally to the 
power spectrum on all scales, i.e. we expect to obtain 
Cn = constant from random point sources. The amplitude 
of this contribution clearly depends on the source counts of 
the population, and can be calculated as 

rScut 

Ct{v) = J ^-S 2 dS„ (23) 

(see e.g. Tegmark & Efstathiou 1996; Scott & White in 
preparation). Then we convert flux units to thermodynamic 
temperature units, as discussed in Section 3.3. Note, how¬ 
ever that here we are interested in the flux rather than flux 
per beam, and so we use the same conversion, without the 
beam solid angle. In equation (23) v is the particular fre¬ 
quency under consideration, S v is the flux, dN/dS^ is the 
differential source count (i.e. N(S„) is the number of sources 
per steradian fainter than S v ), and S cu t is some flux above 
which we feel confident we can identify individual sources. 
All reasonable source counts give convergent Ct at the faint 
end, and generally the precise position of the upper flux cut 
will not be critical. Integrating by parts, and switching to 
the more conventional notation N = N(>S„), the number 
density brighter than some flux threshold, we obtain 
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Flux cut (mJy) 


Figure 9. Qfl at estimates based on the unresolved source count 
model described in the text. The solid line is the normalization 
based on the 5 sources detected in the Hubble Deep Field by 
Hughes et al. (1998), with the 90 per cent confidence region en¬ 
closed by the dashed lines. 


C t (v) = 2 N S v dSv — N(S cut ) Self (24) 

Jo 

In order to estimate what we might expect, we have 
performed this integral using a phenomenological fit to the 
current SCUBA counts. We find that a simple two power-law 
fit of the form 

—©"> ir 

fits all the available data (Smail, Ivison & Blain 1997; Barger 
et al. 1998; Holland et al. 1998a; Hughes et al. 1998, Lilly 
et al. in preparation), and has the same shape as some suc¬ 
cessful models for star formation history (e.g. fig. 11.(b) of 
Blain et al. 1998), with A ~ few x 10 6 Sr -1 , So ~ lOmJy, 
a ~ 0.8 and (3 ~ 1.8. For definiteness we normalise to the 
(Hughes et al. 1998) counts from the Hubble Deep Field 
(HDF), since this seems to be currently the most robust. 
However, the small numbers of sources found in any of the 
deep surveys allows a wide range of possible normalizations. 
A total of 5 sources were detected in the HDF above 3 mJy, 
and so the 90 per cent confidence range for Poisson statis¬ 
tics spans the range 2.0-10.5. We use the central value of 
A = 4.1 x 10 6 Sr -1 to give a specific model for estimating 
the expectations, keeping in mind that a factor of 3 lower or 
2 higher would not be very unlikely. Although our parame¬ 
terized model is a reasonable one, we wish to stress that it is 
just a convenient example; In reality the shape is only con¬ 
strained over a narrow range of fluxes, and of course could 
be quite different at either the bright or faint ends. 

We plot the resulting power spectrum in Fig. 9 as a 
function of the flux cut for bright sources. Explicitly, we 
have converted to the equivalent amplitude of a flat power 
spectrum through the experimental window function (see 


Section 4.2.1), as measured by the expectation value of the 
quadrupole Qns.t- The dashed lines on the plot represent the 
90 per cent CL range from the HDF normalization. We see 
that an equivalent Qa at of higher than around 200 fi.K would 
be expected for even the most drastic cut of 2 mJy. Since 
this is close to the rms in the data-set that we analysed, our 
expectation would be that we could set a limit no lower than 
about 200 /rK. 

For sources below the break, S <C So the integral be¬ 
comes analytic and 


Cc 


Aol / So \ a <-,2 

2 ~~^ VSc© cut ’ 


(26) 


which is just (a/(2 — a))NSc Ut . Hence for a particular nor¬ 
malization of the counts, and for small values of the slope, 
the value of Qfl at will vary approximately as %fa. So anal¬ 
ysis of fluctuations in SCUBA data-sets, such as the one 
presented here, are sensitive to the slope of the counts at 
the faint end. 

The phenomenological model gives a central value of 
267 [i K, even for a flux cut as low as 2 mJy. We note that if 
we were to normalise the counts to effectively higher values, 
such as the Smail et al. (1997) central value, we would pre¬ 
dict fluctuations about 1.3 times higher. Recognising that we 
should have seen evidence of fluctuations in the data if there 
were a large number of faint sources, we can use our limit 
to constrain the faint end slope. To do this we fix the flux 
cut to be 2 mJy, which is very conservative for this data-set, 
and also convenient since it corresponds to the flux cut for 
the HDF counts (Hughes et al. 1998). Then we can fold our 
likelihood distribution for Qfl a t together with the Poisson 
probability distribution for the 5 HDF sources to calculate 
a likelihood function for the faint end slope. We find a 95 per 
cent confidence limit of a < 0.52, which is already shallower 
than some of the models. We stress that this is quite con¬ 
servative in terms of assuming the lowest feasible flux cut 
for our data-set. More stringent constraints could possibly 
be placed by calculating counts and fluctuations for detailed 
models, but this is probably not warranted for the current 
modest data-set. 

The above discussion considered only sources that are 
distributed randomly on the sky. On large scales it is rea¬ 
sonable to assume that the clustering of distant sources 
does have a negligible effect. This is because any three di¬ 
mensional clustering will be washed out when projected 
in two dimensions. However, this is no longer justifiable 
at the smallest angular scales probed, and particularly for 
the sources which may dominate fluctuations detectable by 
SCUBA, namely dusty star-forming galaxies at high red- 
shift. In hierarchical clustering scenarios we expect galaxies 
and clusters to have formed from the build up of smaller 
units at earlier times. The angular scales probed by SCUBA 
lie in the range which is typical for a rich cluster of galaxies 
at high redshift. Hence we might expect to see dusty galaxies 
clustering together as they form a rich cluster, or sub-clumps 
of large galaxies being assembled on group scales, or perhaps 
galaxies flaring up in star formation as they fall into groups 
or are gobbled up by centrally dominating cluster galaxies. 
It is to be hoped that issues such as these may be addressed 
using detailed fluctuation analyses of the deepest SCUBA 
integrations. 
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6 CONCLUSIONS 

The data-set analysed here was obtained using less than an 
hour of integration time, which makes it significantly shorter 
than most other CMB experiments! In order to improve on 
these results, a deeper integration is necessary. The average 
signal-to-noise for the binned data points in the likelihood 
fits is roughly 0.1, although we erred on the side of under¬ 
binning. Most CMB experiments aim for a S/N of unity, 
which is generally considered to be the optimal compromise 
between sky coverage and integration time. This is certainly 
feasible with more ambitious SCUBA integrations, including 
existing data-sets. It should be pointed out however, that 
there may be additional challenges involved with applying 
correlation analyses to data taken in the mapping mode. 

Using the 450 /j,m channel as an atmospheric monitor 
does introduce more noise in the final results, but at levels 
of only a few percent. This indicates that another method 
of sky removal is possible, if there is good reason not to use 
the 850 pm data-set itself. There is also some evidence of 
residual atmospheric gradients across the SCUBA array, but 
they are small, and do not appreciably change our results. 
With the growing amount of blank field SCUBA data, both 
these effects can and should be studied in more detail. 

We have presented here the best upper limit of CMB 
anisotropy on arcsecond scales at 850 /j,m. A careful treat¬ 
ment of the data was performed in order to test for 
non-astronomical signals. Our best estimate for the likeli¬ 
hood function yields an upper limit on the flat-spectrum- 
extrapolated quadrupole of Qnat < 152 /iK. Using a recent 
estimate of the source counts at 2mJy we can convert this 
into an upper limit on the slope of the faint counts, since oth¬ 
erwise we would have detected the fluctuations due to even 
Poisson-distributed sources. We find that if N(> S) oc S~ a 
then a < 0.52 at the 95 per cent confidence limit. 

Larger SCUBA data-sets already exist, and these could 
easily be subjected to similar analysis techniques to those 
discussed here. We fully expect that future data-sets will 
yield detections, since the fluctuations expected in most 
favoured models lie only just below what we could have 
detected here. Investigation of such signals should provide 
independent constraints on the population of sources which 
are just below the detection threshold, and should also yield 
new insight into the clustering properties of the SCUBA- 
bright objects. It is also of course possible that other phys¬ 
ical effects will be important at these wavelengths and an¬ 
gular scales, for example Sunyaev-Zel’dovich effects, or per¬ 
haps currently unconsidered foreground processes within our 
own Galaxy. In any case, fluctuation studies with SCUBA 
are likely to be increasingly important for understanding 
the sub-millimetre and far-infrared sky, and in particular 
for planning future satellite missions such as FIRST and 
Planck. 
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